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Abstract 

In this paper we investigate the effectiveness of Alternating Direction Imphcit (ADI) time 
discretization schemes in the numerical solution of the three-dimensional Heston-HuU-White 
Qh , partial differential equation, which is semidiscretized by applying finite difference schemes on 

' nonuniform spatial grids. We consider the Heston-HuU-White model with arbitrary correla- 

tion factors, with time-dependent mean-reversion levels, with short and long maturities, for 
cases where the Feller condition is satisfied and for cases where it is not. In addition, both 
European-style call options and up-and-out call options are considered. It is shown through 
extensive tests that ADI schemes, with a proper choice of their parameters, perform very well 
in all situations - in terms of stability, accuracy and efficiency. 
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^ : 1 Introduction 



X 



The main aim of this paper is to investigate the effectiveness of Alternating Direction Implicit 
(ADI) time discretization schemes in the numerical solution of three-dimensional time-dependent 
partial differential equations (PDEs) arising in financial option valuation theory. As a prototype 
case we consider the Heston-Hull- White PDE, but our conclusions concerning ADI schemes extend 
to many other related three-dimensional models. 

Consider the asset price process given by the system of stochastic differential equations (SDEs) 

dSr = RrSr CLt + Sr dW^ , 

dVr = Kir] - K) dr + aiVWdW^ , (1.1) 
dRr = a{b{T) - Rr) dr + 02 dW^ . 

The random variables S-n Vr, Rt represent, respectively, the asset price, its variance and the 
interest rate at time r > 0. The parameters k, rj, tJi and a, 02 are given positive real constants 
and h denotes a given deterministic, positive function of time, called the mean-reversion level. The 
VF^, VF^, are Brownian motions under a risk- neutral measure with given correlation factors 
Pill Pi3j P23 G 1] such that the pertinent correlation matrix is positive semideflnite. 

The asset price model (|l.ip can be viewed as an extension of the popular Heston stochastic 
volatility model (Heston (1993)) where the interest rate is not constant but also follows a stochastic 
process, described here by the Hull-White model (Hull & White (1990)). The function b is chosen 
as to match the current term structure of interest rates. The hybrid Heston-Hull-White model 
(jl.ip has recently been studied in Giese (2006), Muskulus, In 't Hout, Bierkens et al (2007), 
Grzelak, Oosterlee & Van Weeren (2009), Grzelak & Oosterlee (2011) and can lead to a more 
accurate valuation of option products that are sensitive to both volatility and interest rates. 
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Let T > be given. If at time re [0, T) the asset price equals s, the variance equals v and the 
interest rate equals r, then for a European-style option with maturity time T and payoff function 
(j) the risk-neutral value is given by 

ifis, V, r, r) = E [e" ''^ (/)(S't, Vt, Rt) | S'r = s, K = = r-j , (1.2) 

where E denotes conditional expectation under the risk-neutral measure. In this paper we consider 
w(s, w, r, t) = ip(s, V, r, T—t). Common arguments in financial mathematics imply that if the option 
value function u is sufficiently smooth then it satisfies the PDE 

du 12 '5^^ 1 2 ^^■^ 1 2^^^ 

cp'ii O^u d^u 

+ Pl2<JlSV—-—- + Pl3(J2sVv—-— + P23Cri^2\/^- 



SOT asar ovor 

du du du 

+ rs— — hK(fi — v)- — h a(6(r — t) — r)— ru (1.3) 

OS ov or 

for s > 0, u > 0, — oo < r < oo and < t < T. Here — oo < r < cxo since the Hull-White 
model yields any, positive or negative, value for the interest rate. We refer to (|1.3p as the Heston- 
Hull-White (HHW) PDE. It forms a time-dependent convection-diffusion-reaction equation on 
an unbounded, three-dimensional spatial domain. The HHW PDE contains three mixed spatial- 
derivative terms, stemming from the correlations between the underlying Brownian motions. Next, 
if w J, then all second-order derivative terms, apart from the d^u/dr'^ term, vanish. This degen- 
eracy feature is already familiar from other financial PDEs, such as the Heston PDE. Finally, we 
note that the coefficient of the du/dr term is time-dependent. 

The HHW PDE is complemented by initial and boundary conditions that are determined by 
the specific option under consideration. The initial condition is given by the payoff function, 

u(s, V, r, 0) = (j){s, V, r). (1-4) 

Boundary conditions will be discussed below. 

The initial-boundary value problem for the HHW PDE does not admit analytic solutions in 
(semi) closed-form in general. An exception concerns European call options if the two correlations 
Pi3 and P23 are equal to zero. Then a direct extension of Heston's (1993) formula is available; it 
is given in the Appendix. 

For the numerical solution of the HHW PDE we consider the well-known and versatile method- 
of-lines approach, see eg, Hundsdorfer & Verwer (2003). Here the PDE is first discretized in the 
spatial variables s, v, r. This leads to a system of stiff ordinary differential equations, the so- 
called semidiscrete system, which is subsequently solved by applying a suitable time discretization 
method. Since the HHW PDE is three-dimensional, the obtained semidiscrete systems are very 
large and also possess a large bandwidth. As a consequence, the selection of the time discretization 
method is critical for its effective numerical solution. To this purpose, we analyze in the present 
paper splitting schemes of the ADI type. 

An outline of the rest of our paper is as follows. In Section[2]we describe the spatial discretiza- 
tion of the HHW PDE. Here finite difference schemes on nonuniform spatial grids are applied. In 
Section [3] we formulate and discuss the four ADI schemes under consideration in this paper: the 
Douglas scheme, the Craig-Sneyd scheme, the modified Craig-Sneyd scheme and the Hundsdorfer- 
Verwer scheme. In Section|3]extensive numerical tests with these ADI schemes are presented. Here 
we investigate in detail the temporal discretization errors. Our tests include arbitrary correlation 
factors, time-dependent mean-reversion levels, cases where the Feller condition is satisfied and 
cases where it is not. In addition, both European call options and up-and-out call options are 
considered. Section [5] gives conclusions and issues for future research. 
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2 Space discretization of the HHW PDE 



In this section we describe the spatial discretization of the HHW PDE. For ease of presentation, 
we consider here European cah options. Thus 4'{s,v,r) = max(0,s — K) with given strike price 
K > Q. The spatial discretization is readily adapted to various exotic options; cf also Section 



2.1 Boundary conditions 

For the semidiscretization, the spatial domain is first restricted to a bounded set [0, S'max] x 
[OjKnax] X [-i?max, -Rmax] with fixcd values ^max, Kiax, -Rmax choscu Sufficiently large. The fol- 
lowing boundary conditions are imposed, 

u(s, u, r, i) ~ whenever s = 0, (2.1a) 

— {s,v,r,t) = 1 whenever s = S'max, (2-lb) 
OS 

u{s,v,r,t) = s whenever w = Knax, (2-lc) 
du 

— (s,v,r,t) ~ whenever r = ii^max- (2-ld) 
or 



Clearly these conditions are of Dirichlet and Neumann type. Condition (|2.1b ) is obvious, (|2.1b ) 
and (|2.1b ) have already been used in the literature for the Heston PDE, and (|2.m ) appears to be 
new. Concerning the latter condition, it is straightforward to prove that under the Black-Scholes 
model the rho of a European call option vanishes for extreme values of the spot interest rate, and 
it is plausible that this holds under the asset price model (11.11) as well. 

At the important, special boundary w = we consider inserting v = into the HHW PDE0 
This is motivated by a theorem of Ekstrom & Tysk (2011) revealing that in the Cox-Ingersoll- 
Ross model, which corresponds to Vr in the SDE (|l.ip . the resulting equation is fulfilled by the 
risk-neutral option value. We note the remarkable fact that this holds irrespective of whether or 
not the Feller condition 2k?7 > crj, well-known from the SDE literature, is satisfied. 



2.2 Spatial grid 

The HHW PDE is semidiscretized on a nonuniform Cartesian spatial grid. The nonuniform grid 
defined in this section is advantageous over a uniform one. This will be illustrated by numerical 
experiments in Sectional 

In the s-direction we consider placing relatively many mesh points throughout a given interval 
[S'loft, bright] C [0, S'max] Containing the strike K. This is natural, firstly, because this is the region of 
interest in applications, and secondly, it alleviates numerical difficulties due to the initial (payoff) 
function </) that has a discontinuous derivative aX, s — K. Let integer m\ > 1 and parameter di > 
and let equidistant points ^min ^ < £.1 < ■ ■ ■ < 6rii ~ 6nax be given with 

^ • 1,-1 / "S'loft 

^min = smh 



di 

Srie-ht ~ S'loft 



^max = 6nt + siuh 



dl 

1 / 5'max ^ Slight 



dl 



Note that ^min < < ^int < Cmax- The mesh = sq < si < . . . < s,„j = S'max is then defined 
through the transformation 

Si = vi^i) (0 < i < mi) 



^We are grateful to Peter Forsyth for a stimulating discussion on this issue. 



3 



S !E- xxmmixx x 



200 



400 



600 



800 



1000 



1200 



1400 



:keo<x XXX 



10 



r 



X X xmsxx X X 



-1 -0.8 -0.6 -0.4 -0.2 0.2 0.4 0.6 0.8 1 

Figure 1: Sample meshes for s, v, r with mi = m2 = ms = 20 and K = 100, T = 1, c = 0.1. 



where 

' 5left + dl Sinh(0 (^min < C < 0), 

= < ^left + dl? (0<e<eint), 

+ dl Sinh(^ - ^int) (6nt < € < €max)- 

This mesh for s is uniform inside the interval [S'left) •S'right] and it is nonuniform outside. The 
parameter di controls the fraction of points Si that lie inside. Put = ~ Co- It is readily seen 
that the above mesh is smooth, in the sense that there exist real constants Co, Ci, C2 such that 
the mesh widths Asj = Sj — Sj-i satisfy 

Co < Asi < Ci A^ and \Asi+i - Asi| < C2 {A^f (uniformly in i, mi). 

For the v- and r-directions we define nonuniform meshes of the same type as considered in, 
eg, Tavella & Randall (2000) and In 't Hout & Foulon (2010). Let integers m2, ma > 1 and 



parameters c, ^2, ^3 > and let equidistant points ?7o < ??i < 
be given by 

Vj = i • A?7 (0 < i < m2) 



< ?7m2 and Co < Ci < • • • < Crr 



with 

and 
with 



Then meshes 
defined by 



AC 



Ar]= — sinh ^(Kiax/(i2), 

•~12 

- c)/d3) + A; • AC (0 < fc < ma) 

[sinh"^((i?i„ax - c)/d3) - sinh"\(-i?max - c)/d3)] . 



Cfe = sinh-^((-A 
1 



ms 

Vo < Vi < ... < Vm2 



To < ri < . . . < 



Rrr 



are 



Vj = d2 smh{T]j) (0 < j < m2) and rfc = c + ds sinh(Cfc) (0 < A; < ms). 

It is easily verified that the meshes for v and r defined above arc also smooth. The parameters 
d2 and ds control, respectively, the fraction of points Vj that lie near v = and the fraction of 
points Vk that lie near a given interest rate level r = c. Here c is chosen depending on the specific 
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mean-reversion function b. For the w-mesh, besides the fact that the region u w is of practical 
importance, it is natural to place relatively many mesh points there for numerical reasons, as 
the HHW PDE is convection-dominated in the u-direction for w w and the initial function is 
nonsmooth. 

In this paper we set S'max = 14ii', Vmax — 10, Rmax = 1- This renders the error induced by the 
restriction of the spatial domain of the HHW PDE to be negligible in our experiments. Based on 
numerical tests, the parameters of the grid have been taken equal to di ~ K/20, d2 = V^max/500, 
da = i?max/400 and, with r = i, 

^left = max{i, e"'''^}i^ , S'l-ight = K. 

A further investigation into possibly better parameter values than above may be interesting, but 
this is out of the scope of the present paper. Figure [I] displays sample meshes for the three spatial 
directions if mi — m2 = = 20 and K — 100, T — 1, c — 0.1. It is clear that the mesh points 
in the s-, v- and r-directions are concentrated, respectively, near s ^ K , v ~ and r = c. 

2.3 Finite difference discretization 

Let / : M R be any given function, let {xi}i^z be any given increasing sequence of mesh points, 
and Axi — Xi — Xi-i for all i. To approximate the first and second derivatives of /, we employ 
the following well-known FD formulas: 

f'{xi) « a-2f{xi-2) + a-if{xi-i) + aof{xi), (2.2a) 
f'{x,) « /?-i/(x,_i)+/?o/(xO+/3i/(a;.+i), (2.2b) 

f'{x^) ~ l0I{x^)+llf{x,+l)+^2f{x^+2). (2.2c) 

f"[xi) « 5_i/(x,_i) + (5o/(x,) + <5i/(x,+i) (2.2d) 



with 



•^-2 — A2;._i(Aa;._i+Aa;i) ' 0.-1 — Aaj^-iAx. ' ~ Axi{Axi-i+Axi)'' 

3 _ -Axj+i a _ Aaj+i-Azj n _ Axj 

1 ~ Axi{Axi+Axi + i) ' "0 — Aa:iAa;i+i ' "1 ~ Axi-^i{Axi+Axi^i) ' 

—2Axi^i — Axi^2 Axi^i+Axi^2 —^^i + i 

^0 — Aa;i+i(A2:i+i-|-Aa:i+2) ' ~ Axi-i-iAxi+2 ' ^'^ ~ 2^772(3^1+7+3^172)' 

-2 r. _ 2 



r _ 2 r _ —2 r 



ArEi(A3:i+Aa;i-|.i) ' AsiAxi + i ' ^ Aa;i+i(Aa;i+Aa;i+i) ' 

Note that (j2.2b[) and (|2.2d|) are central schemes whereas (|2.2al) and (|2.2cp are backward and 
forward schemes, respectively. If / : — R is any given function of two variables (x, y), then we 
approximate the mixed derivative fxyixi,yj) at any point {xi,yj) by successive application of the 
scheme (j2.2b[) in the x- and y-directions. This is equivalent to a FD formula based on a 9-point 
stencil centered about {xi,yj). The FD schemes under consideration all possess a second-order 
truncation error on smooth meshes whenever / is sufficiently often continuously differentiable. 

The actual FD discretization of the initial-boundary value problem for the HHW PDE is 
performed as follows. In view of the Dirichlet conditions (|2.1b .) and (j2.1b ). the relevant set of grid 
points is 

g = {{si,Vj,rk) :l<i<TOi,0<j<m2-l,0<fc< 7713}. 

At this grid, each spatial derivative appearing in (|1.3p is replaced by its corresponding central FD 
approximation, except: 
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• In the region v > rj the backward scheme (|2.2al) is apphed for du/dv. This is done to 
alleviate spurious oscillations in the FD solution when ci is small. It is well-known that 
such oscillations notably arise with central schemes if there is strong advection towards a 
Dirichlet boundary. 

• At the boundary s = S'max the derivatives in the s-direction need to be considered. The 
Neumann condition p.lb ) of course yields du/ds and it subsequently implies that d^u/dsdv 
and d^u/dsdr vanish there. Next, d^u/ds^ is approximated at s = Smi = ^max by the 
scheme (j2.2d|) with virtual point + As,n^ > S'max where the value at this point is defined 
by linear extrapolation, using the value at Smi-i and the (given) derivative at Smi- 

• At the boundary u = we consider setting u = in the HHW PDE, see Subsection 12.11 
Here du/dv is approximated using the forward scheme (|2.2cp . We remark that this is done 
independently of whether or not the Feller condition holds. All other derivative terms in the 
ii-direction vanish if v — and are trivially dealt with. 

• At the boundaries r — ±i?niax the Neumann conditions (j2.1t i) are incorporated similarly as 
for s above. 

The FD discretization of the initial-boundary value problem for the HHW PDE leads to an 
initial value problem for a system of stiff ordinary differential equations (ODEs), 

U'{t) = A{t)Uit) + g{t) (0 < t < T), [/(O) = Uq. (2.3) 

Here A(t), for < t < T, is a given real square matrix and g{t) is a given real vector that is deter- 
mined by the boundary conditions. The entries of the solution vector U{t) form approximations 
to the option values u(s, v, r, t) at the spatial grid points (s, v, r) e Q, ordered in a convenient way. 
The vector C/(0) = Uq is directly obtained by evaluation of the initial function at Q. 

We refer to (12. 3|) as the semidiscrete HHW PDE. The size of this system equals M = 
TOim2(m3 -|- 1) and is very large in general. In the experiments in this paper, we shall deal 
with sizes up to approximately one million. 

3 Time discretization: ADI schemes 

Selecting a suitable time discretization scheme for the semidiscrete HHW PDE (|2.3p is the key to 
obtaining an effective full numerical solution method for the HHW initial-boundary value problem. 
Popular standard methods such as the Crank-Nicolson scheme are often not efficient anymore. The 
reason for this lies in the fact that in each new time step very large systems of linear equations 
need to be solved involving the matrix A(t) for one or more new values of t. Due to its large 
bandwidth, this is computationally very demanding. 

For the time discretization of the semidiscrete HHW PDE, we consider in the present paper 
splitting schemes of the ADI type. Here, the matrix A{t) is decomposed into four simpler matrices, 

A{t)=Ao + Ai + A2 + A3it). 

The matrix Aq represents the part of A{t) that stems from the FD discretization of all mixed 
derivative terms in the HHW PDE. Note that Aq is nonzero whenever at least one of the correlation 
factors pi2, pi3, P23 is nonzero. In line with the classical ADI idea, the matrices Ai, A2, A^lt) 
represent the parts of A{t) that stem from the FD discretization of all spatial derivatives in the s-, 
V- and r-directions, respectively. The ru term in (|1.3p is distributed evenly over Ai, A2, A3{t). We 
decompose 5(t) = go+gi+g2+g3{t) analogously to A(i). The matrices Ai, A2, A3{t) are essentially 
tridiagonal, pentadiagonal and tridiagonal, respectively. Note that the time-dependency of A{t) 
is only passed on to the matrix A3{t), ie, the matrices Ag, Ai, A2 are time- independent. 

Let 9 > he a given real parameter and At = T/N with integer > 1. Set i„ = n/S.t and 
^9n = .93 (in) — 93 (in-i) = .9 {tn) — 9 (^n-i)- Wc study four ADI schemes which all generate. 
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in a one-step manner, successive approximations [/„ to the solution vectors U{tn) of 
n= l,2,...,Af. 



for 



Douglas (Do) scheme: 

' = Un-l + At {A (tn-l) Un-l + 9 (tn-l)), 

Yj = Yj-i + e At A, {¥, - Un-l) (j = 1, 2), 

Yz^Y2+6At (^3 {tn) Yi - A3 {tn-l) Un-l + Ag„ 
Un = ilj. 



(3.1) 



Craig-Sneyd (CS) scheme: 



Ya 


= Un 


-1 + At {A {tn-l) Un-l + g {t, 


.-1)), 








i+eAtAj{Y,~Un-i) (j = 


-1,2), 






= Y2 


+ 9 At {A3 {tn)Y3^A3 {tn-l) 


Un-l 4 


- Ag„) 


% 


= Ya 


+ ^AtAo{Y3-Un-l), 






Y, 


= Y,- 


i+eAtAj{Yj~Un-i) (j = 


-1,2), 






= Y2 


+ 61 At {A3 {tn)Y3~A3 {tn-l) 


Un-l 4 


- Ag„) 


Un 


= ?3 









(3.2) 



Modified Craig-Sneyd (MC'S) scheme: 



Yo 


^Un 


-1 + At {A {tn-l) Un-l 


+ 5(t«-i)), 




Y, 


= Y,- 


i+eAtA,{Y,-Un-i 


) (J = l,2), 




Y3 


= Y2 


+ 6 At {A3 {tn)Y3~A3 


{tn-l) Un-l + 


A5n), 


% 


^Yo 


+ eAtAo{Y3~Un-l), 






Yo 


= Yo 


+ {^-0) At {A{tn)Y3 


~A{tn-l)Un 


-1 + Agn 


Y, 


= Y,- 


i+eAtA,{Y,-Un-i 


) (J = l,2), 




Y3 


= ?2 


+ eAt {A3 {tn)Y3~A3 


{tn-l) Un-l + 


A5„), 


Un 


-?3 









(3.3) 



Hundsdorfer-Verwer (HV) scheme: 



Yo 


= Un 


_1 + At {A {tn-l) Un-l +9 {tn-l)), 


Y, 


= Y,- 


l+9AtAj{Yj-Un-l) (j=l,2), 


Y3 


= Y2 


^6 At {A3 {tn) Y3 - A3 {tn-l) Un-l + Agn), 


Yo 


= Yo 


+ I At {A {tn) Y3~A {tn-l) Un-l + Agn) , 


Y, 


= Y,- 


i+dAtA,{Yj-Y3) (i-1,2), 


% 


= ?2 


^9 At A3 {tn){Y3-Y3), 


Un 







(3.4) 
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The CS, MCS, HV schemes can be viewed as different extensions to the Do scheme. The CS 
and MCS schemes are equivalent if (and only if) = 

It is readily observed that in the four ADI schemes the Aq part, representing all mixed deriva- 
tives in the HHW PDE, is always treated in an explicit fashion. The first papers to propose this 
kind of adaptation of the classical ADI schemes to PDEs with mixed derivative terms are, to our 
knowledge, McKee & Mitcheh (1970) and Craig & Sneyd (1988). 

Following the classical ADI approach, the Ai, A2, A3,{t) parts are treated in an implicit fashion. 
In every step of each scheme, systems of linear equations need to be solved, successively involving 
the matrices {I — 6 At Aj) for j = 1,2 and {I ~ 6 At A3(t„)), where / is the identity matrix. As all 
these matrices have a fixed, small bandwidth (of at most five) this can be done efficiently by LU 
factorization. Note that for j = 1,2 the pertinent matrices are further independent of the step 
index n, and hence, their LU factorizations can be computed once, beforehand, and then used in 
all time steps. 

By Taylor expansion one obtains (after some elaborate calculations) the classical order of 
consistency of each ADI scheme, ie, the order of consistency in the nonstiff sense. For any given 9, 
the order of the Do scheme is just one if ^0 is nonzero. This low order is due to the fact that 
the Aq part is treated in a simple, explicit Euler fashion. The CS scheme has order two provided 

= ^. The MCS and HV schemes are of order two for any given 6. With the latter schemes, the 
parameter 9 can thus be chosen to meet additional requirements. 

A detailed discussion, with ample references to the literature, concerning the above four ADI 
schemes has been given in In 't Hout & Welfert (2007, 2009). The Do and CS schemes are already 
often applied to PDEs in finance, see eg, Andersen & Andreasen (2000), Lipton (2001), Randall 
(2002) and Andersen & Piterbarg (2010). More recently, the MCS and HV schemes have gained 
interest, see eg. In 't Hout (2007), Dang, Christara, Jackson & Lakhany (2010), In 't Hout & 
Foulon (2010), Haentjens & In 't Hout (2010), EgloflF (2011) and Itkin & Carr (2011). 

For an effective application of numerical schemes, stability is imperative. The stability of ADI 
schemes in the case of PDEs possessing mixed derivative terms has been analyzed by a number of 
authors in the literature. This stability analysis has been performed in the von Neumann (Fourier) 
framework. Here one considers application to the semidiscretized convection-diffusion equation 

du 

— = c • V-it + V • (DWu) 
ot 

on a rectangular domain, with constant real vector c and constant, positive semidefinite real matrix 
D ~ {dij), with periodic boundary condition, on a uniform spatial grid, and one studies stability 
in the ^2-iiorm. Note that the presence of mixed derivative terms corresponds to the matrix D 
being nondiagonal. A desirable property is unconditional stability, ie, without any restriction on 
the time step Ai > 0. 

The most comprehensive stability results for the Do, CS, MCS and HV schemes in the literature 
up to now, relevant to PDEs with mixed derivative terms, are given in In 't Hout & Welfert (2007, 
2009), In 't Hout & Mishra (2010, 2011). We review the main conclusions from loc cit pertinent 
to two and three spatial dimensions. Here stability is always understood in the von Neumann 
sense and unconditional. To formulate some of the results, we consider for 7 e [0, 1] the following 
condition on _D, 

\dij \ < 7 ^diidjj for aU i ^ j. (3.5) 

The quantity 7 can be viewed as a measure for the relative size of the mixed derivative coefficients. 
Because D is positive semidefinite, the condition p.5|) is always fulfilled with 7 = 1. But in actual 
applications, in particular the HHW PDE, one usually has more information, namely 7 < 1. 

For two-dimensional convection-diffusion equations with mixed derivative term, the Do and 
CS schemes are both stable whenever 9 > ^. If there is no convection (c — 0), then the MCS 
and HV schemes are stable whenever 9 > ^ and 9 > 1 — (~ 0.293), respectively. For the 

MCS scheme, stability has been proved for general two-dimensional equations, with convection, if 

1 < 6* < 1. Next, based on strong numerical evidence, stability of the MCS scheme for the special 
value = \ was conjectured under the mild, additional condition that p.Sp holds with 7 < 0.96. 
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For the HV scheme, stabihty for general two-dimensional equations has been conjectured for all 
6* > 5 + (« 0.789). We note that the latter bound stems from Lanser, Blom & Verwer (2001), 
who proved it to be necessary and sufficient for stability in the case of two-dimensional equations 
without mixed derivatives. 

For three-dimensional problems, positive results on the stability of the ADI schemes have been 
derived for pure diffusion equations with mixed derivative terms. In this case, it has been shown 
that the Do, CS, MCS and HV schemes are stable whenever 6* > |, 6* > i, 6* > max{i, ^(27 + 1)} 
and 61 > 1(2 - VS) (« 0.402), respectively!! 

At this moment sufficient conditions on 9 for stability of the ADI schemes pertinent to general 
three-dimensional convection-diffusion equations with mixed derivative terms are lacking in the 
literature. Accordingly, we select the parameters 9, in the subsequent experiments, on the basis 
of the present results, reviewed above. 

In practical applications it turns out that a smaller value 9 often leads to a smaller error 
constant. In view of this, we choose 9 as small as possible under the requirement of (unconditional) 
stability. 



4 Numerical experiments 

In this section we present extensive numerical tests with the four ADI schemes (13. ip . p.2p . p.3p . 
p.4p in the application to the semidiscrete HHW PDE described in Section [2l This yields im- 
portant insight in their actual stability and convergence behavior and mutual performance. We 
consider the HHW model with arbitrary (nonzero) correlation factors, with time-dependent mean- 
reversion levels, for cases where the Feller condition is satisfied and for cases where it is not. In 
addition, we deal with European call options as well as up-and-out call options. 
For the diffusion matrix of the HHW PDE, 

D{s,v)^-\ pi2(Jisv afv P2z(yxa2\Jv , 

V P\Z(y2S^/v /5230-lO-2v^ ^2 / 

it is easily verified that the condition p.Sp holds with 7 = max{|pi2|, Ipial, |p23|}- Based on the 
stability and accuracy results discussed in Section [3] we select, for this value 7, 

• the Do scheme (|XT|) with 6* = | 

• the CS scheme with 9 = \ 

• the MCS scheme with 9 = max{i, ^(27 -I- 1)} 

• the HV scheme ([Q]) with 9 ^\^\\pi. 

The Do scheme has classical order one and the CS, MCS, HV schemes all possess classical order 
two. Note that the MCS scheme has 6* < ^ (« 0.462). 

For the ADI schemes under consideration we shall study in this section the global temporal 
discretization error, defined by 

e{At;Tni,m2,ni3) = max{\Ui{T)~UN,i\- {si,Vj,rk) e V}, (4.1) 

where T = NAt with integer > 1 and U (T) denotes the exact solution vector to the semidiscrete 
HHW PDE (|2.3p at time T. The index I = l{i,j, k) corresponds to the spatial grid point (s^, Vj,r}S) 
and P is a natural region of interest, to be specified below. 

If P13 = P23 = 0, then a semi closed-form analytic formula for European call option values is 
known, see the Appendix. We shall employ this formula to validate the FD discretization of the 

^The result for the Do scheme is new; its proof will be included in a forthcoming paper. 
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Case A 


Case B 


Case C 


Case D 


Case E 


Case F 


K 


3 


0.6067 


2.5 


0.5 


0.3 


1 


V 


0.12 


0.0707 


0.06 


0.04 


0.04 


0.09 




0.04 


0.2928 


0.5 


1 


0.9 


1 


a 


0.2 


0.05 


0.15 


0.08 


0.16 


0.22 


Cl 


0.05 


0.055 


0.101 


0.103 


0.055 


0.074 


C2 


0.01 


0.005 


0.001 


0.003 


0.025 


0.014 


C3 


1 


4 


2.3 


1 


1.6 


2.1 


a-2 


0.03 


0.06 


0.1 


0.09 


0.03 


0.07 


Pl2 


0.6 


-0.7571 


-0.1 


-0.9 


-0.5 


-0.3 


Pl3 


0.2 (0) 


0.6 (0) 


-0.3 (0) 


0.6 (0) 


0.2 (0) 


-0.5 (0) 


P23 


0.4 (0) 


-0.2 (0) 


0.2 (0) 


-0.7 (0) 


0.1 (0) 


-0.2 (0) 


T 


1 


3 


0.25 


10 


15 


5 


K 


100 


100 


100 


100 


100 


100 



Table 1: Parameters for the Heston-Hull- White model. 



HHW PDE from Section [5] and to study the global spatial discretization error in this case, defined 
by 

e(mi,m2,m3) = max{ |u(si,Wj,rfe,r) - J7;(T)| : (si, Wj, rfc) e P }. (4.2) 

The temporal and spatial discretization errors are both measured in the maximum norm, which 
is highly relevant to financial applications. In order to compute (14.11) and (|4.2I) for a given spatial 
grid, we use a sufficiently accurate reference value for U{T), obtained by applying the MCS scheme 
to (1231) with N = 20000 and N = 200 time steps, respectively. 

For efficiency of the spatial discretization it turns out that one can place relatively less grid 
points in the v- and r-directions than in the s-direction. Accordingly, we choose in the following 
the numbers of grid points in the three spatial directions as mi = 2to, TO2 = TO3 = m with 
integer m. Note that the size of the semidiscrete HHW system equals Af = 2m^{m + I). 

We are interested in mean-reversion levels b that are time-dependent. As an example, we 
choose 

6(t) = Cl - C2e-^-^" (r > 0) (4.3) 

with positive constants ci, C2, C3 and ci > 02- This choice for b is somewhat arbitrary, but the 
conclusions obtained below on the numerical schemes are the same for other (more realistic) time- 
dependent mean-reversion levels. For the mesh in the r-direction, defined in Subsection 12.21 we 
take c — Cl. 

4.1 European call options 

Our first experiments concern European call option values in the six cases of parameter sets for 
the HHW model hsted in Tabled! 

The cases A, B, C can be viewed as an extension of three test cases for the Heston model 
previously used in In 't Hout & Foulon (2010). The values k, ry, cti, pi2, T stem from Bloomberg 
(2005), Schoutens, Simons & Tistaert (2004) and Winkler, Apel & Wystup (2002), respectively. 
Here the Feller condition always holds. 

The cases D, E, F form an extension of the three cases for the Heston model presented by 
Andersen (2008). They are proposed in loc cit as challenging test cases for practical applications. 
Notably, the Feller condition is not fulfilled. Also, the maturity times are large. 

In all six cases, the values a, ci, C2, C3, (72 pertinent to the Hull- White model as well as the two 
correlations pia, P23 are chosen in an arbitrary, realistic way. Here the corresponding correlation 
matrices are always positive definite. 
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Figure 2: Spatial discretization errors e(2m, m, m) vs l/m for European call options in the six 
cases of Table [T] with pi3 = p23 = for m = 10, 15, ... , 75. 
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Figure 3: Spatial discretization errors with uniform 2m x m x m grid for European call options 
in the six cases of Tabled] with pi3 = p23 = for m ~ 10, 15, ... , 75. 
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We first consider the FD discretization and study the spatial discretization errors defined by 
(|4.2p . with region of interest 

I?=(iX,fi^)x (0,l)x (0,i). 

As mentioned above, it is possible to compute these whenever pi3 = /323 = 0. Figure [2] displays the 
errors e(2m, m, m) vs 1/m in the six pertinent cases of Table [T] for m = 10, 15, ... , 75. Note that 
m — 75 means M — 855000 spatial grid points, which was the practical (memory) limit on our 
laptop computer. Figure [2] clearly shows that in each case the spatial discretization errors decrease 
as m increases. To determine the numerical orders of convergence, straight lines have been fitted 
to the results. In the cases A, B, C, D, F the obtained orders of convergence are all equal to two 
approximately. Only in case E a slightly lower order was obtained, namely 1.6. As an indication 
of the sizes of the spatial discretization errors in a relative sense, we mention that these always 
lie between 0.2% and 1.2% when m = 50 and between 0.1% and 0.6% when m = 75 (here only 
option values are considered that are greater than 1). In view of the foregoing, wc conclude that 
the FD discretization defined in Section [5] performs satisfactory in all six cases. It is interesting 
to briefly compare the spatial errors to those obtained with a uniform grid and the same number 
of grid points. Figure [3] shows spatial discretization errors analogously to Figure [21 but then for 
uniform 2m x m x to grids. Clearly, in most cases the spatial errors for the nonuniform grid are 
substantially smaller, often by an order of magnitude, than those for the corresponding uniform 
grid. Further, it is clear that for a uniform grid the behavior of the spatial error as a function 
of the number of grid points is erratic, which is undesirable. Also, the nonuniform grid yields 
more points in the region in (s, v, r)-space where one wishes to obtain option prices. We therefore 
conclude that the nonuniform grid defined in Section [2] is preferable over a uniform grid. 

We next consider the performance of the four ADI schemes in the application to the semidiscrete 
HHW PDE for European call options in the six cases of Table [1] with all correlations nonzero. 
Figure m displays the temporal discretization errors e(At;2m,m,m) for a sequence of step sizes 
with 10^3 < At < 10° when m = 50. 

A first main observation from Figure |4] is that for all four ADI schemes the temporal discretiza- 
tion errors are bounded from above by a moderate value and decay monotonically as At decreases. 
Additional experiments indicate that this is true for any value to; see a further discussion below. 
This suggests an unconditionally stable behavior of the schemes, which is a new and nontrivial 
result. It does not directly follow for example from the von Neumann stability analysis presented 
in Section [31 We note that this result holds in all six cases, independently of whether or not the 
Feller condition is fulfilled. 

A next observation is that the CS scheme exhibits an undesirable feature in the cases A, B, C 
with temporal errors that are very large for moderate At, compared to what may be expected on 
the basis of its asymptotic error behavior (ie, for small At). To a much lesser extent, this is also 
observed with the HV and MCS schemes. Additional experiments reveal that the relatively large 
temporal errors occur at spatial grid points near the strike K. It is already known in the literature 
that the nonsmoothness of the initial function at the strike yields high-frequency errors which are 
not always sufficiently damped by numerical schemes, notably the Crank-Nicolson scheme and 
the Do and CS schemes with 9 = ^. A popular remedy for this situation is to first apply, at 
t = 0, two implicit Euler steps with step size At/2, and then to proceed onwards from t = At 
with the scheme under consideration, cf Rannacher (1984). However, in our present application of 
the three-dimensional HHW PDE this damping procedure is computationally intensive. We shall 
consider an alternative in the next subsection. 

A further analysis of the results in Figured indicates that in each case the temporal discretiza- 
tion errors for the Do scheme are bounded from above by CAt and for the MCS, HV schemes by 
C(At)^ (whenever At > 0) with constants C depending on the scheme and the case. This clearly 
agrees with the respective orders of consistency of the schemes. Moreover, experiments with both 
smaller and larger values of to suggest that the constants C are only weakly dependent on the 
number of spatial grid points M, ie, the error bounds are valid in a stiff sense, which is a desirable 
property. This result is also nontrivial, as the order of consistency is a priori only relevant to fixed, 
nonstiff systems of ODEs. For the CS scheme, we find that the temporal errors can be bounded 
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Figure 4: Temporal discretization errors e (Ai; 100, 50, 50) vs At for European call options in the 
SIX cases of Tabled ADI schemes: Do with 6* = § (diamond), CS with 6* = i (dark circle), MCS 
with 9 = max{i, ^(27+1)} (light circle) and HV with 6* = i + i%/3 (square). 
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in each case by C(Ai)^ with a constant C independent of stiffness if damping is apphed. Actual 
numerical experiments for ADI schemes combined with damping will be presented in the next 
subsection. 

Our implementation of the ADI finite difference discretization has been done in Matlab, where 
all matrices have been defined as sparse. For the CS, MCS, HV schemes the cpu-time per time 
step was about 0.10, 0.18, 0.90, 1.5 cpu-seconds for m = 25, 30, 50, 60, respectively, on one Intel 
Core Duo T7250 2.00 GHz processor with 4 GB memory; for the Do scheme these times are about 
halved. Here all correlations were nonzero and the mean reversion level was time-dependent. It 
readily follows that the cpu-times are indeed almost directly proportional to the number of spatial 
grid points M ^ 2m^. 

4.2 Up-and-out call options 

As an important and particularly challenging type of exotic options we consider here European- 
style up-and-out call options. The FD discretization described in Section [5] is adapted with few 
modifications. Let barrier ^max ~: B > K he given. Then the boundary conditions (12.1b ). (|2.1b ) 
are replaced by 

u{s,v,r,t) = whenever s = B, (4.4a) 
du 

— {s,v,r,t) = whenever = Kiax- (4.4b) 
ov 

The condition (|4.4b ) has been suggested by various authors in the literature. Note that all 
boundary conditions are now homogeneous, and g(t) = 0. The relevant set of spatial grid points 
is 

Q = {{si, Vj, ^k) : 1 < * < ™i ~ 1 , < J < m2 , < fc < ma}. 

The only significant change we make to the FD discretization of Section [5] is to replace, in the 
s-direction, the central advection scheme (|2.2b[) by the backward scheme (|2.2ap if r < and by the 
forward scheme (|2.2cP if r > 0. This upwind approach alleviates spurious oscillations in the FD 
solution that are obtained with the central advection scheme. It is already useful for up-and-out 
call options in the one-dimensional Black-Scholes model. The pricing of up-and-out call options is 
numerically more challenging than of vanilla options, due to the boundary layer that is introduced 
at the barrier. 

Figure [S] displays the numerically obtained up-and-out call option values in the six cases of 
Table[l]for barrier B = 120 and (sampled) spot interest rates r « 0.02 on the (s, w)-domain [0, B] x 
[0, 1). Here the FD discretization has been applied with m = 50 and for the time discretization 
the HV scheme is used with At = 10~^. 

We study in detail the performance of the four ADI schemes. Similar to the case of European 
call options. Figure IHl shows the temporal discretization errors e(At;2TO,m,m) in the case of 
up-and-out call options for a sequence of step sizes 10~^ < At < lO'^ when m — 50. As a first 
observation, it is clear from Figure[6]that the unfavorable feature of relatively large temporal errors 
for moderate step sizes is more pronounced compared to the case of vanilla options, especially for 
the CS scheme, cf Subsection 14.11 We attribute this to the additional discontinuity of the payoff 
function at the barrier s = B. We therefore consider application of a damping procedure at t = 0. 
Instead of performing two substeps at t = with step size At/2 by the implicit Euler scheme, 
which forms a common approach, we employ here the Do scheme, with parameter value 9 = \. 
This is computationally more attractive when dealing with multidimensional PDEs. Figure [7] 
shows the temporal discretization errors in the case of up-and-out call options when two initial 
substeps with the Do scheme and 9 — 1 are applied. Clearly, the behavior of the temporal error as 
a function of the step size has become regular and, in most cases, at only a limited loss of accuracy 
for small At (an exception being the CS scheme in case E). Hence, the present damping procedure 
performs satisfactory. Applying two substeps of the implicit Euler scheme for the damping would 
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Case A 



Case B 




Figure 5: European up-and-out call option values in all cases of Table [T] with barrier B ~ 120. 
Spot interest rates in the cases A, B, C, D, E, F are, respectively, r = 0.025, 0.022, 0.025, 0.027, 
0.022, 0.017. Note: scales on vertical axes vary. 
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Case A Case B 




Figure 6: Temporal discretization errors e (At; 100, 50, 50) vs At for up-and-out call options in all 
cases of Table [T] for barrier B — 120. ADI schemes: Do with — ^ (diamond), CS with & = 
(dark circle), MCS with 6 = max{i, ^(27+1)} (light circle), and HV with 6 = 5 + ^\/3 (square). 
No initial damping. 
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Figure 7: Temporal discretization errors e{At; 100, 50, 50) vs At for up-and-out call options in all 
cases of Table [T] for barrier B — 120. ADI schemes: Do with d — j (diamond), CS with — ^ 
(dark circle), MCS with 6 = max{i, ^(27+1)} (light circle), and HV with 6 =\ + \^/i (square). 
Two initial damping substeps using the Do scheme with 6 = 1. 
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yield similar or somewhat smaller temporal errors than those in Figure [T] However, we find that 
this comes at a much higher computational cost, also when iterative solvers, like BiCGSTAB, are 
applied. We thus infer that damping with the Do scheme is more efficient. 

As a main positive conclusion, the numerical results for all ADI schemes are consistent with 
an unconditionally stable behavior: the temporal discretization errors are bounded from above by 
a moderate value and decay monotonically as At decreases, which is obtained for any value m 
tested. A closer inspection of the results displayed in Figure [7] yields that the temporal errors 
behave for sufficiently small At as C(At)P with p = 1.0 for the Do scheme and 1.6 < p < 2.0 for 
the CS, MCS, HV schemes, with constants C. Experiments with different values of m reveal that 
both p and C are only weakly dependent on the number of spatial grid points M, indicating that 
the error behavior is valid in a stiff, hence favorable, sense. Note further that the difference in 
performance between the Do scheme and the CS, MCS, HV schemes is often less striking than in 
the case of vanilla options, but for the latter three schemes combined with damping still a higher 
order and increased accuracy is obtained. 

5 Conclusions and future research 

In this paper we studied ADI schemes in the numerical solution of the three-dimensional HHW 
PDE: the Do scheme, the CS scheme, the MCS scheme and the HV scheme, each with a well chosen 
parameter 9. Extensive experiments have been conducted for six cases of parameter sets for the 
HHW model, including correlations that are all nonzero, time-dependent mean-reversion levels, 
and short and long maturities. In three cases the Feller condition is not fulfilled. We consid- 
ered both European call options and up-and-out call options. Our tests have shown that all ADI 
schemes perform very well in terms of stability, accuracy and efficiency. In particular they always 
reveal an unconditionally stable behavior. Next, the Do scheme always has a stiff order of con- 
vergence equal to one. The CS, MCS, HV schemes show a stiff order of convergence equal to two 
for European call options and, when combined with damping, between 1.6 and 2.0 for up-and-out 
call options. 

Based on the numerical experiments and the theoretical stability results, we find that the MCS 
scheme with 9 = max{i, ^(27 -I- 1)} and the HV scheme with 9 — ^ + are preferable. Here 
7 = maxij \pij\. Also the CS scheme with 6* = i is a good candidate. For the latter scheme, a 
damping procedure at t = is always recommended. Damping can be done efficiently, in an ADI 
fashion, by applying the Do scheme with 9 = 1. 

The Do, CS, MCS and HV schemes are expected to perform well and possess similar favorable 
properties as obtained in this paper in the numerical solution of many other three-dimensional 
PDEs and for other exotic options. Also, the ADI schemes can directly be applied, with high 
efficiency, when any other FD discretization is employed, as the matrices Aj{t) {1 < j < 3) always 
have a small bandwidth. We shall investigate other applications in future research. At the same 
time, a further theoretical stability analysis of the ADI schemes will be carried out. 

Appendix 

Here wc give the semi closed-form analytic formula for European call option values tp{s,v,r,T) 
under the HHW model (|1.1|) with pis = p23 = as derived in Muskulus, In 't Hout, Bierkens et 
al (2007). The notation is adapted to our present situation. We put p = pi2. 

The solution presented in loc cit is of a form similar to the Black-Scholes formula, 

(p{e^ , V, r, r) = e^Pi{x, v, r, r) — KB{r, t)P2{x, v, r, r). 

Here B{r,T) denotes the value at time r of a zero-coupon bond that pays 1 at maturity, given 
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that at time r the short rate equals r. For this, it is well-known that 



B{r,T) 
c{r,T) 
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The Pi, P2 can be viewed as probabilities and are retrieved from characteristic functions /i, /2 
by inversion: 
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The above valuation formula is easily seen to constitute a proper extension of Heston's (1993) 
formula, by taking 6(r) = rg and (T2 = 0. It can be approximated to any accuracy, by a direct 
adaptation of numerical integration techniques already well studied in the literature for Heston's 
formula. We note that the additional integrals involving the function h can be exactly determined 
in our particular case of 1 



Acknowledgements 

The authors gratefully acknowledge Peter Forsyth for a stimulating discussion, convincing them 
of the proper boundary condition for the HHW PDF at w = and pointing them to the work 
by Fkstrom & Tysk. The authors also thank Jan Van Casteren for a valuable, unpublished note 
on the validity of the HHW PDF. Furthermore, they are indebted to Sven Foulon for providing 
an implementation of Heston's formula for Furopean call options, which we extended to our case. 
This work has been supported financially by the Research Foundation - Flanders, FWO contract 
no. G.0125.08. 



20 



References 



L. Andersen, Simple and efficient simulation of the Heston stochastic volatility model, J. Comp. 
Finan. 11 (2008) 1-42. 

L. Andersen & J. Andrcasen, Jump- diffusion processes: volatility smile fitting and numerical 
methods for option pricing, Rev. Deriv. Research 4 (2000) 231-262. 

L. B. G. Andersen & V. V. Piterbarg, Interest Rate Modeling, Volume I: Foundations and 
Vanilla Models, 1st ed., Atlantic Financial Press, 2010. 

4] Bloomberg Quant. Finan. Dcvel. Group, Barrier options pricing under the Heston m,odel, 2005. 

51 I. J. D. Craig & A. D. Sneyd, An alternating- direction implicit scheme for parabolic equations 
with mixed derivatives, Comp. Math. Appl. 16 (1988) 341-350. 

D. M. Dang, C. C. Christara, K. R. Jackson & A. Lakhany, A PDE pricing framework for 
cross- currency interest rate derivatives, Proc. 10th Int. Conf. Comp. Sc. (ICCS), Proc. Comp. 
Sc. 1 (2010) 2371-2380. 

D. Egloff, GPUs in financial computing part III: ADI solvers on GPUs with application to 

stochastic volatility, Wilmott mag., March 2011, 51-53. 

E. Ekstrom & J. Tysk, Boundary conditions for the single-factor term structure equation, Ann. 
Appl. Prob. 21 (2011) 332-350. 

A. Giese, On the pricing of auto-callable equity structures in the presence of stochastic volatility 
and stochastic interest rates, presentation MathFinance Workshop, Frankfurt (2006). 
Available at www.mathfinance.com/workshop/2006/papers/giese/slides.pdf 

10] L. A. Grzclak & C. W. Oosterlce, On the Heston model with stochastic interest rates, SIAM 
J. Finan. Math. 2 (2011) 255 286. 

11] L. A. Grzelak, C. W. Oosterlee & S. van Weeren, Extension of stochastic volatility equity 
models with the Hull-White interest rate process, published online in Quant. Finan. (2009), 
doi:10.1080/14697680903170809. 

12] T. Hacntjcns & K. J. in 't Hout, ADI finite difference discretization of the Heston Hull White 
PDE, In: Numerical Analysis and Applied Mathematics, eds. T. E. Simos et al, AIP Conf. 
Proc. 1281 (2010) 1995-1999. 

13] S. L. Heston, A closed-form solution for options with stochastic volatility with applications to 
bond and currency options. Rev. Finan. Stud. 6 (1993) 327-343. 

14] K. J. in 't Hout, ADI schemes in the numerical solution of the Heston PDE, In: Numerical 
Analysis and Applied Mathematics, eds. T. E. Simos et al, AIP Conf. Proc. 936 (2007) 10-14. 

15] K. J. in 't Hout & S. Foulon, ADI finite difference schemes for option pricing in the Heston 
model with correlation, Int. J. Numer. Anal. Mod. 7 (2010) 303-320. 

16] K. J. in 't Hout & C. Mishra, A stability result for the Modified Craig-Sneyd scheme applied 
to 2D and 3D pure diffusion equations. In: Numerical Analysis and Applied Mathematics, 
eds. T. E. Simos et al, AIP Conf. Proc. 1281 (2010) 2029-2032. 

17] K. J. in 't Hout & C. Mishra, Stability of the modified Craig-Sneyd scheme for two-dimensional 
convection-diffusion equations with mixed derivative term. Math. Comp. Simul. 81 (2011) 2540- 
2548. 



18] K. J. in 't Hout & B. D. Welfert, Stability of ADI schemes applied to convection- aiTtusion 
equations with mixed derivative terms, Appl. Numer. Math. 57 (2007) 19-35. 



21 



[19] K. J. in 't Hout & B. D. Welfert, Unconditional stability of second- order ADI schemes applied 
to multi- dimensional diffusion equations with mixed derivative terms, Appl. Numer. Math. 59 
(2009) 677-692. 

[20] J. Hull & A. White, Pricing interest-rate- derivative securities, Rev. Finan. Stud. 3 (1990) 
573-592. 

[21] W. Hundsdorfer & J. G. Verwer, Numerical Solution of Time-Dependent Advection-Diffusion- 
Reaction Equations, Springer, Berlin, 2003. 

[22] A. Itkin & P. Carr, Jumps without tears: a new splitting technology for barrier options, Int. 
J. Numer. Anal. Mod. 8 (2011) 667-704. 

[23] D. Lanser, J. G. Blom & J. G. Verwer, Time integration of the shallow water equations in 
spherical geometry, J. Comp. Phys. 171 (2001) 373 393. 

[24] A. Lipton, Mathematical Methods for Foreign Exchange, World Scientific, Singapore, 2001. 

[25] S. McKee & A. R. Mitchell, Alternating direction methods for parabolic equations in two space 
dimensions with a mixed derivative. Computer J. 13 (1970) 81-86. 

[26] M. Muskulus, K. J. in 't Hout, J. Bierkens, A. P. C. van der Ploeg, J. in 't Panhuis, F. Fang, 
B. Janssens & C. W. Oosterlee, The ING problem: a problem from, the financial industry, Proc. 
58th European Study Group Mathematics with Industry, eds. R. H. Bisseling et al, Utrecht 
(2007) 91-115. 

[27] C. Randall, PDE Techniques for Pricing Derivatives with Exotic Path Dependencies or Exotic 
Processes, Lecture notes. Workshop GANdiensten, Amsterdam, 2002. 

[28] R. Rannacher, Finite element solution of diffusion problems with irregular data, Numer. Math. 
43 (1984) 309-327. 

[29] W. Schoutens, E. Simons & J. Tistaert, A perfect calibration! Now what?, Wilmott mag., 
March 2004, 66-78. 

[30] D. Tavella & C. Randall, Pricing Financial Instruments, Wiley, New York, 2000. 

[31] G. Winkler, T. Apel & U. Wystup, Valuation of options in Heston's stochastic volatility model 
using finite element methods, in: Foreign Exchange Risk, eds. J. Hakala & U. Wystup, Risk 
Books, London (2002) 283-303. 



22 



